!-------------NEW subroutines (V. Stepanenko, May 2007): rad_par2, eps1, epsil, sin_sun

      SUBROUTINE rad_par2(taixiy,qa,pres,tclds,atm_rad,net_solar)
!-----Computes solar and atmospheric radiation according to simple
!-----parameterizations      
      use alloc
      implicit none
      real(8), external:: eps1,ndaym,epsil
      real(8) sigma,coef_l,tclds,atm_rad,net_solar,sol_const,
     & taixiy,qa,hume,pres,days,delta,tha,sinh0,phi1,h0,
     & coef_s,a_ang,d_ang,c_ang,a_brent,b_brent	 
      REAL(8), EXTERNAL:: sin_sun
      integer(4) i
      logical firstcall
      data firstcall /.true./
      SAVE
      
      if (firstcall) then
        sigma = 5.67d-8
        coef_l = 0.22
        sol_const = 1376.
        coef_s = 1. - 0.4393
!       tau_rg = 2.74 !3.5 !0.4 
 
        a_ang = 0.82
        d_ang = 0.25
        c_ang = 0.95/100.

        a_brent = 0.526
        b_brent = 0.065/10.
      endif

      hume = qa*pres/0.622
      select case(rat_par)
        case(1)
          atm_rad = taixiy**4*sigma*eps1(taixiy,hume)
        case(2)
          atm_rad=sigma*taixiy**4*(a_ang-d_ang*10.**(-c_ang*hume))
        case(3)
          atm_rad=sigma*taixiy**4*(a_brent+b_brent*dsqrt(hume))
      end select
      atm_rad = atm_rad*(1+coef_l*tclds*tclds)

      sinh0 = sin_sun(imonth,iday,ihour,iminu,iseco,xlatit)
      h0 = dasin(sinh0)
      net_solar = sol_const*(1-coef_s*tclds)*sinh0/
     & (1+epsil(h0)*tau_rg/sinh0)
      if (sinh0 < 0) net_solar = 0
 
      if (firstcall) firstcall=.false.
      RETURN
      END SUBROUTINE rad_par2
